Astronomy & Astrophysics manuscript no. paper 


© ESO 2008 


February 5, 2008 





Collisional dust avalanches in debris discs 



Anna Grigorieva 1 , P. Artymowicz 1,2 , and Ph. Thebault 1,3 

\Q ■ 

O , Stockholm Observatory, SCFAB, SE-10691 Stockholm, Sweden 

^SJ 2 University of Toronto at Scarborough, 1265 Military Trail, Toronto, Ontario, MIC 1A4, Canada 
Observatoire de Paris, Section de Meudon, F-92195 Meudon Principal Cedex, France 



<D ■ 

{J*} ' Received / accepted 

<n : 



ABSTRACT 



> 

in : 
o ■ 

On 1 

o : 

O ■ 

^ : 

Oh 1 
I 1 

o 1 

a: 



We quantitatively investigate how collisional avalanches may develop in debris discs as the result of the initial breakup of a 
planetesimal or comet-like object, triggering a collisional chain reaction due to outward escaping small dust grains. We use a 
specifically developed numerical code that follows both the spatial distribution of the dust grains and the evolution of their 
size-frequency distribution due to collisions. We investigate how strongly avalanche propagation depends on different parameters 
(e.g., amount of dust released in the initial breakup, collisional properties of dust grains, and their distribution in the disc). Our 
simulations show that avalanches evolve on timescales of ~ 1000 years, propagating outwards following a spiral-like pattern, and 
that their amplitude exponentially depends on the number density of dust grains in the system. We estimate a probability for 
witnessing an avalanche event as a function of disc densities, for a gas-free case around an A-type star, and find that features 
created by avalanche propagation can lead to observable asymmetries for dusty systems with a /3 Pictoris-like dust content or 
higher. Characteristic observable features include: (i) a brightness asymmetry of the two sides for a disc viewed edge-on, and (ii) 
a one-armed open spiral or a lumpy structure in the case of face-on orientation. A possible system in which avalanche-induced 
structures might have been observed is the edge-on seen debris disc around HD 32297, which displays a strong luminosity 
difference between its two sides. 
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1. Introduction 



Direc t imaging: of ci r cumst e llar discs ( e.g.. iHeap et al 
2000; IClampin et all l2003t IlTtI 12004 | S chneider et al 



2005) have provided resolved disc morphologies for several 



systems (e.g., /3Pic, HD 141569A, HD 100546, HD 32297) 
and have shown that dust distribution is not always 
smooth and axisymmetric. Warps, spirals, and other 
types of asymmetri es are commonly observed (e.g., 
iKalas fe Jewittlfl995l for the /3Pic system). These mor- 
phological features can provide hints on important ongo- 
ing processes in the discs and improve our understanding 
of the evolution of circumstellar discs and of planetary 
formation. 

The usual explanation proposed for most of these 
asymmetries is the perturbing influence of an embedded 
planet. As an example, the warp in the /3Pic disc has 
been interpr eted as induced by a jovian planet on an in - 
clined orbit l|Mouillet et al.lll997tlAugereau et al.ll200lj) . 
Likewise, for annulus-like discs with sharp inner or outer 
edges, the most commonly proposed explanation is trun- 
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cation or gap o pening due to planets or bou nd stellar 
companions (e.g.. lAueereau fc Papaloizoul2004h . although 
alternative mechanisms such as gas drag on dust grains 
within a gas disc of limited exte nt have also been proposed 
l|Takeuchi fe Artvmowic3l200l[) . For spiral structures, au- 
thors have also been specul ating on gravitational instabil- 
ities ( Fukagawa et al. 2004 1), as well as on a bound stellar 
companion l|Augereau fc Papaloizou1l2004l) . 

The catastrophic breakup of one single large object re- 
leasing a substantial amount of dust fragments could be 
an alternative explana tion for some observed asymmetries. 
IWvatt fe DentJ (j2002) have examined how such collision- 
ally produced bright dust clumps could be observed in 
Fomalhaut's debris disc. Likewise , such clumps have been 
proposed by iTelesco et all l|2005l) as a possible explana- 
tion for mid-infrared brightness asymmetries in the cen- 
tral f3 Pictoris disc, but only based on preliminary order of 
magnitude estimate s . Mor e recently, the detailed study of 
iKenvon fc Bromley! ((2005) investigated the possibility of 
detecting catastrophic two-body collisions in debris discs 
and found that such a detection would require the breakup 
of 100-1000 km objects. The common point between these 
different studies is that they focus on global luminosity 
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changes due to the debris cloud directly produced by the 
shattering events themselves. 

In the present paper, we re-examine the consequences 
of isolated shattering impacts from a different perspective, 
i.e., by considering the collisional evolution of the pro- 
duced dust cloud after its release by the shattering event. 
The main goal here is to s tudy one pos s ibly v ery efficient 
process, first proposed bv lArtvmowic d (|l097f) . but never 
quantitatively studied so far, i.e., the so-called collisional 
avalanche mechanism. The basic principle of this process 
is simple. After a localized disruptive event, such as the 
collisional breakup of a large cometary or planetesimal- 
like object, a fraction of the dust then produced is driven 
out by radiation pressure on highly eccentric or even un- 
bound orbits. These grains moving away from the star 
with significant radial velocities can breakup or micro- 
crater other particles farther out in the disc, creating in 
turn even more small particles propagating outwards and 
colliding with other grains. Should this collisional chain 
reaction be efficient enough, then a significant increase in 
the number of dust grains could be achieved. In this case, 
the consequences of a single shattering event, in terms of 
induced dust production, could strongly exceed that of 
the sole initially released dust population. The outward 
propagation of the dusty grains could then induce observ- 
able asymmetric features in the disc, even if the initially 
released dust cloud is undetectable. 

The goal of this work is to perform the first quan- 
titative study of the avalanche process and investigate 
the morphology of avalanches in debris discs, under the 
assumption that dust dynamics is not controlled by gas 
( Lagrange et alJ EoOO"). For this purpose we have created 
a numerical code, described in Sect. that enables us 
to simulate the coupled evolution of dynamics and size- 
frequency distribution of dusty grains. The results of our 
simulations, which explore the effect of several parameters 
(total mass and radial distribution of dust in the disc, mass 
and size distribution of the planetesimal debris, physical 
properties of the grains and the prescription for collisional 
outcome for grain-grain collisions) are presented in Sect.0] 
In Sect. we examine under which conditions avalanche- 
induced features might become observable. We end with 
a discussion of the probability of witnessing an avalanche 
(Sect-EJ) and finally a summary (Sect.EJ). 



2. Simplified theory of dust avalanches 

A dust avalanche is a chain reaction of outflowing debris 
impacting disc particles and creating even more debris ac- 
celerated outwards by the star's radiation pressure. The 
basic principle of this mechanism can be illustrated by a 
set of analytical equations. We present here a simplified 
theory of avalanches based on the order-of-magnitude ap- 
proach of lArtvmowiczl l|l997f ). firstly for its pedagogical 
virtues, but also because it can serve as a reference that 
facilitates the understanding of the main results derived 
from our extensive numerical exploration. 



Let us assume that N particles of size s gr (radius) 
move through a cloud of dust grains of size s at a rela- 
tively high velocity Let us further assume that each colli- 
sion produces a constant number Np of such debris, which 
are quickly accelerated to velocities leading to further de- 
structive collisions. To derive the total number of debris 
produced by the avalanche, we define the optical depth as 



dr = n(s)a(s)dl, 



(1) 



where n(s) is the number density of dust particles of size 
s in the system, a(s) rs tt(s + s gr ) 2 is the cross-Sect, for 
collisional interaction between grains, and dl is the length 
measured along the grain path. The number of debris pro- 
duced in the interval dr is then 



dN = NN P dr. 



(2) 



Integration over the whole path of the grains gives the 
total number of debris produced by the avalanche 



iVtot = N exp(Npr) 



(3) 



where Nq is the number of outflowing grains initially re- 
leased. 

In a disc, r can be approximated by the optical thick- 
ness in the disc midplane, 



J J ns 2 dn(s)dR, 



(4) 



where R is the radial cylindrical coordinate. We replace 
Np by its average value (Np) to emphasize the fact that 
in reality Np depends on the details of each collision. 
Equation then takes the form 



N tot ~ N exp({Nf))T]\). 



(5) 



This equation gives an estimate of an avalanche efficiency 
in a disc through the total number of grains A^ot it pro- 
duces. However, one should keep in mind that the rele- 
vance of this set of equations is limited to global, order 
of magnitude estimates. Furthermore, these equations do 
not give any insight into the temporal development and 
spatial structure of a given avalanche. For these crucial 
issues, numerical modeling is clearly required. 

3. The model 

The number of dust grains in a circumstellar disc is far 
too large to follow every grain individually during the cal- 
culation; some kind of statistical approach must therefore 
be used. Models of dust disc evolution developed to date 
fall into two main categories. On the one hand, "parti- 
cle in a box" models divide the dust grains into statis- 
tical bins according to their size and enable us to com- 
pute the evolution of the size distr ibution within a give n 
spatially homogeneous region ('e.g.. lThebault et "al1l20n.'^l . 
While it is possible to mimic a spatially inhomogeneous 
system by integrating a set of coupled particle-in-a-box 
models, this can become unwieldy in the absence of strong 
simplifying symmetries. iKenvon fc Bromlevl 1 2004^ use a 



A. Grigorieva et al.: Collisional avalanches in dusty discs 



3 



multiannulus code for example, but their model is one- 
dimensional in space. On the other hand, direct N-body 
simulations (treating the dust as test particles in the po- 
tential of a 2 or multi-body system) are used to accurately 
follow the spatial evolution of dynamical st ructures such 
as planet induced gaps or res onances (e.g., IffiiaS [2003; 
lAugereau fc Panaloizoull2004|) . In this case, however, the 
sizes of the dust grains are either not taken into account 
or assumed to be equal. 

For the present problem, however, we need to follow 
both the spatial distribution of the grains and their size 
distribution with reasonable accuracy. To do this, we de- 
veloped a new code in which all grains with similar param- 
eters (size, chemical composition, spatial coordinates, and 
velocity) are represented by a single superparticle (here- 
after SP). We follow the dynamical evolution of these SPs 
and compute the collisional destruction and production of 
grains as SPs pass through each other. We represent newly 
created grains as new SPs. The maximum number of SPs 
our code can handle is about one million. 



3.1. Superparticles 

A detailed description of our SP modeling is given in the 
appendix. Here we briefly outline its main characteristics: 
A SP is described by the position and velocity of its cen- 
ter of mass (which coincides with its geometrical center), 
by its size, shape, and internal density profile, and by 
the number of dust grains it contains. For the present 
work all SPs are treated as cylinders and their geometri- 
cal centers are constrained to lie in the midplane of the 
disc. The cylinders have constant radii r sp and variable 
heights h sp (R), where R is the distance from the star (see 
Appendix EJ . All grains inside a given SP are assumed 
to have the same physical properties. We assume that all 
grains in our simulation are spheres with identical den- 
sities, chemical compositions, and porosities. The grains 
(and thus the SPs) are distributed into mass bins sep- 
arated by a factor 2 logarithmic mass increment (i.e., a 
factor of 1.26 in size). 

The trajectory of a SP corresponds to the trajectory 
of a test particle (with dynamical properties identical to 
the SP's grains) located at the SP's center of mass (see 
Sect. I3.2JI . SPs can overlap and freely pass through each 
other. In this event, collisional interactions between their 
respective grain populations is considered. This process is 
treated as a passage of two clouds of grains through each 
other (see Appendix lA.3|l . It results in the loss, by destruc- 
tive collisions, of a fraction of the initial grain populations 
and the production of smaller collisional fragments. These 
newly produced debris are placed into newly created SPs 
in accordance with the grain sizes and velocities. In the 
current version of the code, the centers of all SPs move 
in the same plane and the dust distribution is symmetric 
with respect to this midplane. However, the SPs represen- 
tation method could in principle be used to model systems 
with vertical deformations (e.g., warps). 



The size of a SP is fairly large (r sp = 5AU). This 
puts unavoidable constraints on the spatial resolution of 
our simulations and prevents us from modeling processes 
occurring on scales smaller than the SP radius. It would, 
for example, be difficult to model fine resonant structures 
induced by disc-planet interaction. Moreover, the current 
version of the method with a constant value of the SP 
radius is not applicable to collisional evolution in the inner 
regions (^ 20 AU) of debris discs. Although this limitation 
could be overcome by introducing a dependence of the size 
of a SP on the distance to the star (e.g., r sp ex R), we have 
not implemented it in the current version of the code, since 
our main goal here is to model collisional avalanches that 
propagate outwards, inducing observationally significant 
features in the outer (<; 100 AU) regions of the disc. 

The grains inside a SP do not have explicit vertical ve- 
locity components. To check the validity of this assump- 
tion, we have performed test runs, for which an artificial 
vertical velocity dispersion term was added to the pla- 
nar velocity, which showed no significant departure from 
the in-plane velocities case. Note that a vertical velocity 
component is, however, indirectly taken into account by 
the fact that SP heights increase with distance from the 
star (see Appendix IA.3JI . accounting for the geometrical 
dilution of grain spatial densities. 

3.2. SP trajectories 

As has been mentioned earlier, the trajectory of a SP is 
identical to the trajectory of a test particle (with mass, 
size, and chemical composition identical to those of the 
SP's grains) located at the SP's center of mass. Test par- 
ticles move in the gravitational field of a star under the 
influence of the stellar radiation force. The equation of 
motion reads: 



GMr 



rad 



'PR, 



(6) 



where m and r are the mass and position of the test 
particle, G is the gravitational constant, M is the star 
mass, and i^ ra d and i*pR are the radiation pressure and 
Poynting- Robertson drag, respectively. In our simulations 
we can neglect the Poynting-Robertson drag, since it acts 
on a timescale much longer than the time intervals con- 
sidered here. 

The radiation pressure force is expressed as a function 
of the gravitational force through the radiation pressure 
coefficient, j3, as 



■ rad 



— /5-fgrav P 



GMm 



(7) 



The parameter (3 is a function of the stellar luminosity, 
grain size, and op tical properties of the grain material 
feurns et al1ll979tl. We us e the Mie theory code devel- 
oped bv lArtvmowic 1 l|l98Sj) to calculate j3 (Fig. [I). 

A 7 th — 8 th order Runge-Kutta method is used for in- 
tegrating test particles trajectories. Although in the sim- 
ulations presented in this paper the dynamics of the SPs 
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Fig. 1. Ratio of the radiation pressure force to the grav- 
itational force vs. grain size for different grain materials 
and porosities, P, calculated for an A5V ((3 Pictoris-like) 
star. The thick lines represent silicate Mgo.gsFeo.osSiOs . 
The thin lines are used for water ice. The solid lines are for 
solid grains, the dashed lines are for 50% porous grains, 
and the dash-dotted lines are for 80% porous grains. 

is purely Keplerian, we have decided not to use analyti- 
cal solutions since the Runge-Kutta integrator allows for 
an easy inclusion of any additional gravitational (due to 
planetary or stellar perturbers) or dissipative forces (such 
as PR and gas drag). 

3.3. Initial dusty disc structure 

The SP representation method is used to model the initial 
dusty disc structure. The total number of SPs for each 
size bin is chosen so that, at any given location in the 
disc, there are at least 2-5 overlapping SPs to account for 
different dynamical characteristics of grains of this size at 
this location. Each of these overlapping SPs thus differs 
from the others by its local velocity. To model the initial 
dust distribution in the disc we use ~ 5 x 10 4 SPs (test 
runs with larger number of SPs do not lead to significant 
changes in the results). The number density of dust grains 
at a given location in the disc is calculated as the sum of 
the grain densities of the overlapping SPs. 

The archetypical, and still by far the best known, 
debris disc of /JPictoris is taken as a reference sys- 
tem for the initial dusty disc structure. In the present 
study we do not aim to model this particular system 
and just adopt its global properties for the dust distri- 
bution. Alternative dust distributions are also explored in 
Sect. 14.51 For the dust profile in /3Pictoris, we take the re- 



sults of lAugereau et al who numerically derived 

the dust distribution giving the best fit to the resolved 
scattered light images as well as the long- wavelength pho- 
tometric data, as a reference. We assume here that all 
grains are produced from parent bodies on circular orbits 
following the b e st-fit parent body distribution given in 
Au gereau et alJ Ij2001l) . where most of the bodies are lo- 
cated within an extended annulus between 80 and 120 AU, 
with a depletion in the inner <50AU region and a sharp 
drop of the den sity distribution outside 120 A U (see for ex- 
ample Fig. 1 of lThebault fc Augereaul l2005'). Grains with 
small P have almost the same orbits as their parent bodies 
(the biggest grains), while smaller grains (i.e., with higher 
/3) have more elliptic orbits depending on their /3 value. 
The initial number of grains as a function of their size 
follows a classical single power law size-frequency distri- 
bution 



dn = N 0s {s/s )' POB ds 



(8) 



where the power-law coefficient, pos = 3. 5, corresponds t o 
an idealized collisionally evolved system l|Dohnanvll96fl|) . 
The minimum size s m i n for the disc grains is given by 
the radiation pressure blow-out cutoff and corresponds to 
~ 2 /im for the compact grains considered in the nominal 
case. The maximum size s max for the disc grains is taken 
to be 1 cm. Runs with an order of magnitude higher max- 
imum grain size give very similar results while being more 
computationally demanding. At the same time we cannot 
lower s max since millimeter particles make up a few per- 
cent of the disc's optical thickness and their contribution 
starts to be significant for the avalanche development. 

The vertical structure of the disc is expressed in terms 
of the vertical geometrical optical thickness, t±, per unit 
length, z, as 



— — = G T — exp 

dz w 



(9) 



where C T is a normalizing constant, t±(R) = 
J J TTs 2 dn(R, s)dz is the vertical optical thickness of the 
disc at distance R from the star, n(R, s) is the number 
density of dust grains of size s, w(R) is the disc width, 
and p z = 0.7 is a parameter, that determines the shape of 
the vertical profile. The disc width changes with radius as 



w(R) = 0.055R n 



R 

Rm 



(10) 



where R m = 117 AU and p w = 0.75 for most of the runs. 
(Alternative values of p w have been explored in test runs, 
which have shown that results only weakly depend on it.) 

3.4. Collisional outcomes 

Collisions are the crucial mechanism for the development 
of the avalanche phenomenon. The result of a collision, 
in terms of the size- frequency distribution of the debris, 
depends on several parameters: projectile and target ma- 
terials and structures, sizes, impact velocities, and angle 
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of incidence. Since it is not possible to model every col- 
lision in such detail, we have to adopt a simplified al- 
gorithm. We assume that the impact energy of colliding 
bodies, E co \, is equally shared between them. Laboratory 
experiments show that this is the case when both bodies 
are made from ide ntical material regardless of their sizes 
in_e^al|EMi|)- E co i is: 



E, 



MiM 2 w 2 ( 



col 



rcl 



2(Mi + M 2 ) 



where v le \ and are the relative velocity and the 

masses of the colliding bodies. The relative velocity be- 
tween grains is derived from relative velocities between 
SPs after compensating for the artificial Keplerian shear 
induced by the SP finite radius. 

Collision outcomes are traditionally divided into two 
classes: (1) catastrophic fragmentation, when the largest 
remaining fragment, My, is less than half of the parent 
body mass, M, and (2) cratering, when M\f > 0.5M. The 
energy per unit mass that is needed to get Mif = 0.5M 
is called the threshold specific energy Q* . If the spe- 
cific energy Q = 0.5E CO \/M received by a body is more 
than Q* , then the collision leads to catast rophic breakup, 
whereas cratering occurs if Q < O* llFuiiwara etaO 
Il977t iPetit fe Farinellalll993l iBenz fe Asphaudll999j) . O* 
is a function of size for whic h we adopt a classi- 
cal power law dependenc e (e.g., lRvan fc Meloshl Fl998: 
iHousen fc Holsappleil 999'). The collisional response of the 
small objects considered in the present work falls into 
the so-called strength regime, where the target's internal 
strength is the dominant factor, for which 



](s/s y 



PQ 



(11) 



where Qq corresponds to the value of the threshold en- 
ergy for size sn- In t his reg ime, Q* decreases wi t h size . 
IHousen fc Holsapplel l|l990|) and iRvan fc Meiosbl l|l998|) 
present a wide range of values for pq. We cannot di- 
rectly apply their result because we are dealing with 
much smaller sizes. As pointed out in lTielens et al.l 1 19941) 
"for submicron-sized bodies, cracks play little role, and 
the strength of a material approaches the ultimate yield 
strength of t he material" , which corresponds to Q* = 
2 x 10 s ergs/g l(Tielens et al.ll994l and references therein). 
The threshold energy size-dependence most probably has 
a knee in the micron-submillimeter size range, but since 
we do not have any information about the slope change we 
decided not to introduce two additional unknown parame- 
ters into the code, but simply to adopt a slightly shallower 
slope, pq — 0.2, for our calculations. 

To account for the effect of different incidence angles, 
we correct the value of Q* by a correction factor x cr cor- 
responding to an average over all incidence angles 

Q = Qhcadon/^cr, 

where x cr = 0.327 l|Petit & Farinellal Il993h . For both 
catastrophic shattering and cratering prescription, we 
use the approach and the algorithm presented in 



IPetit fc Farinellal |l993) . However, we refine this model by 
assuming that the fragment mass distribution produced 
follows a broken power-law instead of a single-index one: 



dn = N dm 



(m/m s ) qi if m < m s 
(m/m s ) 92 if m > m s 



(12) 



Such a change of slope between the small and large 
fragments, always corresponding to a flattening of the 
slope in the small particle range, is indeed sup ported by 
experimental results (e.g.. iDavis fc Rvanl H9 901 . For the 
choice of values for the power-law indexes (qi, qz) and the 
transition mass m s for the slope change , we take the exper- 
imental studies of lDavis fc RvarJ l|l99f)|) as a reference and 
explore values within the range of possible values obtained 
by these authors. The minimum fragment size is assumed 
to be 0.1 /mi, unless otherwise explicitly specified. 

In our calculations we do not consider changes in the 
orbital parameters of the colliding bodies (i.e., SP), since 
this effect is not important for the present study. There 
are 2 reasons for this: (i) the lifetime of an avalanche (typ- 
ically ~ 10 3 years) is very short from the point of view of 
the global disc evolution, thus we can neglect any changes 
in the disc dynamics caused by mutual collisions between 
the disc particles (i.e., "field SPs" in our simulations); (ii) 
the dynamics for the majority of the avalanche SPs are 
controlled by the radiation pressure. Their orbital param- 
eters are thus determined mostly by their /3 values and 
only weakly depend on the velocities at which these SPs 
are born (as is verified in Sect. 4.2.1 for the first generation 
of avalanche grains). In this respect, taking the velocity of 
the center of mass of the colliding grains as the initial 
velocity for the produced debris is a good approximation 
within the frame of our simulations. 



3.5. Initial planetesimal breakup 

As previously mentioned, we assume that the initial source 
of the collisional avalanche is the breakup of a large, at 
least kilomctcr-sized object. We do not perform a simu- 
lation of the initial shattering event itself, but implement 
a simple parametric prescription for the dust released in 
the breakup. In most runs, we consider a "nominal" case, 
in which Mq = 10 20 g of dust is released in the 0.1 fim to 
1 cm range at Rq = 20 AU from the star, unless otherwise 
explicitly specified. It should be noted that the released 
dust mass Mq is the only relevant parameter for our sim- 
ulations. In this respect, the exact process leading to the 
initial release is not crucial. However, when it comes to 
estimating the probability for such a dust-release event to 
occur (as will be done in Sect. lti.lfl . one has to consider the 
mass Mpb of the parent body whose shattering produces a 
mass M of dust. The ratio M /Mpb is obviously < 1, but 
strongly depends on several poorly constrained parame- 
ters, mainly related to the physics of the shattering event. 
For an idealized case when the largest fragment produced 
has mass M;/ = 0.5Mpe and smaller fragments follow the 
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Dohnanyi "equilibrium" size distribution (dn oc s~ 3,5 ds), 
one gets M PB — 10 24 g. However, laboratory and numeri- 
cal studies as well as observations of ast eroid families (e.g., 
iDavis fc RvanllT99Cll: iTanga et al.l ll999) all point towards 
smaller Mif and steeper size distributions for highly dis- 
ruptive impacts of large objects, with indexes typically 
in the —3.7 to —4 range for the largest > 0.01M;/ frag- 
ments and closer to —3.5 for the smallest ones. Using for 
ex ample the fragmentat ion prescription for large objects 
of iThebault et alJ (2003), we determine that for a typical 
shattering at lkm/s, AIpb — 10 21 g= lOA/o, which cor- 
responds to an object of radius ~ 40 km. We shall thus 
assume a nominal M^/Mp-Q ratio of 0.1 for the discussion 
in Sect. 16.11 For the size spectrum of the dust particles 
released in the 0.1 /im to 1 cm range, we assume a single 
power law (Eq. [5] with p s o = 3.5) for our nominal case. 
The dependence of the results on M , p s o, and other pa- 
rameters related to the planetesimal debris is explored in 
Sect.O 

4. Results 

For the sake of the readability of the results, it is conve- 
nient to divide the system into two populations: 1) the 
avalanche particles, representing all bodies initially re- 
leased by the planetesimal breakup plus all grains later 
created by collisions between the avalanche particles and 
the disc material, and 2) the field particles, i.e., the pop- 
ulation of grains in the disc unaffected by the avalanche 
mechanism. To quantify the magnitude of an avalanche 
we introduce the area amplification factor, F, which is 
the ratio of the total cross-sectional area of the avalanche 
grains, within 500 AU from the star, to the initial cross- 
sectional area of planetesimal debris released. The maxi- 
mum value F max reached by the amplification factor while 
the avalanche is propagating is used to measure the ampli- 
tude of a given avalanche and to compare avalanches ob- 
tained for different initial conditions. Time is expressed in 
orbital periods at 20 AU (~ 70 yr), unless otherwise explic- 
itly specified. Tabled summarizes the set of initial param- 
eters chosen for our "nominal" case. All free parameters 
of the simulations are then explored in separate runs. 



Table 1. Main model parameters for the nominal case. 



Grains: 



material 
porosity 
grain density 



Disc ("field" population): 

minimum size 

maximum size 

radial distribution 

optical thickness along radius 

in the midplane 

disc extension 

Initial planetesimal debris: 
minimum size 
maximum size 
size distribution (see Eq. |8J 
initial mass of dust released 
distance from the star for 
the planetesimal breakup 
initial velocity of 
the center of the mass 

Collisional prescription: 
threshold energy, so = 1 cm 
power-law index (Eq. 111! 
size distribution of debris (see 

minimum size 

power-law indexes for m < m s 
power-law indexes for m > m s 
position of the slope change 



Mg .95Feo.o5Si0 3 
compact grains (P=0) 
3.5g.cm -3 



s min = 2 ^m 
Smax = 1 cm 
Augcrcau e^a^ (20n 



T|| = 0.022 
[20, 500] AU 

Smin.pl = 0.1 flttl 
Smax.pl = 1 Cm 

Po.pi = -3.5 
Mo = 10 20 g 

Ro = 20 AU 

«0 = l.l«kep 

QS = lO^rg.g- 1 
PQ = -0.2 
Ea. lT2jl : 

Smin,col = 0.1 fJ,W. 

qi = 1.5 
q 2 = 1.83 
m s = Mif/3 




5 10 15 20 

Time (in orbital periods at 20 AU) 



25 



4.1. Nominal case (NC) 

Figure |21 shows the temporal evolution of an avalanche, 
for the nominal case, in terms of the vertical opti- 
cal thickness, T±. av , of the avalanche particles (t±^ v = 
JJ irs 2 dn av (s)dz). As expected, the first stages correspond 
to a fast development and multiplication of the avalanche 
grains. In this early expansion phase the surface density is 
dominated by the smallest high-/3 (<; 0.5) particles, which 
contribute to ~ 85% of rj_,av The maximum value of the 
amplification factor is F max = 210 and is reached at t ~ 5 
(~ 350 yrs, see Fig. - After that, the loss of small grains 
on unbound orbits dominates over the collisional produc- 
tion of new dust particles, and the avalanche begins to 
fade. In these later stages, the total cross sectional area of 



Fig. 3. Time evolution of the cross-sectional area amplifi- 
cation factor (the ratio of the total cross-sectional area of 
the avalanche grains within 500 AU to its initial value at 
t = 0). Initial increase is due to dust production by out- 
flowing planetesimal debris colliding with the disc mate- 
rial. When the grain removal (due to star radiation pres- 
sure) rate exceeds the grain production, the value of F 
drops (see text for more details). 



the avalanche grains (within 500 AU) is increasingly domi- 
nated by the bigger grains on bound orbits. It is important 
to point out that the timescale for the avalanche propa- 
gation is short in comparison with orbital periods in the 
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Fig. 2. Nominal case. Color-coded maps (log-scale) of the vertical optical thickness of avalanche grains, Tj_ iav , at 
different stages of the avalanche evolution (t=0.6, 5, 10, 40 orbital periods at 20 AU). The planetesimal debris are 
released at t = at 20 AU from the star. Field particles are not included in the plots. The position of the star is marked 
by the white cross. 
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Fig. 4. Ratio of the total number of grains produced by 
the avalanche by the time t, N to t , to the initial number of 
released planetesimal debris, N . 

outer part of the disc (e.g., only ~ 1/5 of the orbital period 
at 200 AU). 

The amplification achieved by the avalanche mecha- 
nism is impressive, i.e., an increase in grain cross-sectional 
surface density by two orders of magnitude compared 
to the particles initially released by the planetesimal 
breakup (Fig. [3). However, absolute values of T|_ iav are 
still very small compared to those of the field particles, 
with Tj^av/Tj^fieid never exceeding 10~ 2 (see Sect.|S]for a 
more detailed discussion of this crucial parameter). 

To compare the results of our simulation with the sim- 
plified theory of Sect.El we plot the ratio of the total num- 
ber of grains Ntot produced by the avalanche until time 
t to the initial number Nq of released planetesimal debris 
(Fig. 01 . As can be clearly seen, N to t/N quickly reaches a 
plateau, and we take N tot /N ~ 200 at t = 15 as a refer- 
ence value. Plugging values for the average number of par- 
ticles produced by each grain-grain collision, (Np) ~ 150, 
and T|| = 0.022 into Eq.0 we get A r tot,theory/A r o ~ 30, i.e., 
a factor of ~ 7 difference with the result of our simulation. 



This is mainly due to the fact that T|| underestimates the 
real value of t, firstly because the real path of a grain is 
curved rather than parallel to a disc radius, and secondly 
because in T|| the size of the avalanche grains s gr is ne- 
glected. From our simulations, we were able to estimate 
the discrepancy between r and tm to be roughly of a factor 
of 1.6. We thus get 



No 



3 1.6T||(iV3> 



200, 



(13) 



a value close to the numerical results. 

It is also interesting to link N to t/No to the amplifica- 
tion factor parameter F max . By definition, 



F(t) 



Js 2 dN in (s,t) = N in (t){s(t) 2 ) 
Js 2 dN Q (s) N (s 2 ) ' 



(14) 



where N in (s,t) is the total number of avalanche grains 
inside 500 AU at time t, and (s(t) 2 ) and (sg) are the aver- 
aged cross-sectional areas of the avalanche grains at time 
t and of the initially released planetesimal debris, respec- 
tively. Since the avalanche's dust production comes mainly 
from a rather short peak of activity (see Figs. and 0J, 
N to t is close to Ni n (t*), where is the time at which 
Ni n reaches its maximum. A numerical check showed that 
indeed N in (t*) rj 2/3A^ tot , so that 



F„ 



Ntot (s(t*) 2 ) _ 2 (s(U) 2 ) i. a {JV „ )T1 



Nn 



(4) 



(4) 



(15) 



The validity of this relation is easily numerically verified, 
with 2/3A to t/A x (s(i*) 2 )/(s^) ~ 200, a value that is 
indeed relatively close to the -Fmax value obtained in the 
simulation. 
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4.2. Dependence of F m 
debris parameters 



on the initial planetesimal 



4.2.1. Initial mass & velocity of the planetesimal dust 
cloud 

The initial mass Mq of dust released has been explored 
as a free parameter. The simulations show that the max- 
imum amplification factor, Fmax, does not vary with Mq, 
at least in the 10 12 g-10 21 g range, a result that is in good 
agreement with Eqs. El an d El Likewise, F max does not 
change much when varying the initial speed Vq of the cen- 
ter of mass of the planetesimal dust cloud. There is only 
a 20 % increase of F max when vq is increased from Uk ep 
to 1.41«kcp- This weak dependence on the initial velocity 
of the debris confirms the fact that avalanches are driven 
mostly by the smallest particles, which are accelerated to 
high speeds weakly correlated to the initial release veloc- 
ity. 
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Fig. 5. Maximum amplification factor as a function of the 
minimum size assumed for the initial planetesimal debris. 
The power-law index for the size distribution is equal to its 
nominal value po,pi — 3.5. For smaller po,pi the variation 
with s m in,c is less pronounced. 



4.2.2. Size distribution of planetesimal debris 

For the nominal case we choose s m in,pi = 0.1 /urn. This 
value is compatible with the low er limit for the size o f 
the interplanetary dust particles l|Fraundorf et alJll982|) . 
It is also in good agreement with the size distribution de- 
duced from studies of cometary comas that show that 
the sm allest particles are abo u t 0.08 — 0.28 (tm in di- 
ameter l)Kolokolova et alJl200l|) . iMcDonnell et all l)l99l|) 
observed smaller grains in comas, but their contribution 
to the total dust population remained marginal. Even if 
grains smaller than 0.1 /Ltm are produced abundantly in 
the planetesimal breakup, they are not expected to con- 
tribute significantly to the avalanche process since they 
are in the size range where (3 decreases for smaller grains 
(see Fig. nj. As a consequence, they have lower outgoing 
velocities which, together with their smaller masses, lead 
to a marginal contribution in terms of impacting kinetic 
energy. We thus believe 0.1 /im to be a reliable minimum 
value for s m in,pi and explored s m i njP i values in the 0.1 //m 
to 1 /xm range, the latt e r valu e being the one considered 
bv iKenvon fc Bromlevl l)2005l) . Although the impacting 
kinetic energy per grain is increasing in the 0.1 /im to 1 /im 
size range (leading to an increase of the N to t/No ratio), 
-Fmax oc N to t/N x (s(i*) 2 }/ langles^) decreases with in- 
creasing s m in iP i (Fig. 03) because of the decreasing value of 
the (s(£*) 2 )/(so) factor. 

Test runs have also been performed to check the s ma x,pi 
dependence. This exploration has shown that results do 
not depend on this parameter for values higher than 1 cm. 
This is due to the fact that grains bigger than this size have 
very small (3 values, as well as a low total cross-sectional 
area, which do not allow them to significantly contribute 
to the avalanche propagation. 

The dependence of -Fmax on the power-law index for the 
initial planetesimal debris size distribution, po.pi, is shown 
in Fig. El It can be noted that increasing p , P i above 3.5 
(value for the nominal case) does not lead to a significant 
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Fig. 6. Dependence of the maximum amplification factor 
on the power-law index of the initial size-frequency distri- 
bution of the planetesimal debris (Eq.[SJ), s m i n . P i = 0.1 /jm. 



increase of F max , whereas less steep power laws lead to a 
significant decrease in F max . 



4.2.3. Position of the planetesimal breakup 

We perform a set of runs in which the position of the plan- 
etesimal breakup, Ro, is varied between 20 and 100 AU 
(Fig. E}, but all the other parameters remain identical 
to the nominal case. The maximum amplification factor 
decreases with increasing Ro for two reasons: (i) the to- 
tal amount of disc material through which the outflowing 
grains propagate is higher when the grains are released 
close to the star; (ii) the unbound grains (/3 > 0.5) have 
time to reach higher radial velocities if they are released 
closer to the star (see Fig.|HJl, which leads to more violent 
collisions and hence higher dust production per collision. 
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Fig. 7. Dependence of the maximum amplification factor Fig. 9. Dependence of the maximum amplification factor 



on the location of the primary planetesimal breakup. 
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on the minimum size for debris produced in avalanche col- 
lisions (Eq.IHwith qi = 1.5,(72 = 1.83, m s = A/i f /3). 
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Fig. 8. Radial velocities vs. distance from the star for 
grains with different (3 values released by parent bodies on 
circular orbits for a (3 Pic-like star. The release distances 
are 20 AU and 70 AU. 

4.3. Dependence of -F max on the prescription for 
collisional outcome 



Fig. 10. Dependence of the maximum amplification factor 
on the value q\ in Ea. 1121 for collisionally produced grains 
(q 2 = 1.83). 

prevents them from taking up most of the cross-sectional 
area of the avalanche grains. 



4.3.1. Minimum size of avalanche produced debris 4.3.2. Size distribution of the debris grains 

In the nominal case we assume that the minimum size, Numerical exploration of the position of the slope change 

Smin,coi, for the debris produced by collisions and the ™s (with s min = 0.1 ^m) and of the m > m s power-law 

minimum size of the initial planetesimal debris, s min ,pi, index <?2 shows that the resulting amplification factor only 

are both equal to 0.1 fim. We have seen in Sect. weakly depends on these parameters. Changing m s /M u 

that we do not expect significant changes when s mi „, p i < from 1 to 10 leads to changes in F max by only a factor - 2. 

0.1 /mi. However, the situation is slightly different for On the contrary, variation in the gi index can significantly 

avalanche grains, since these grains are continuously pro- affect F max , especially for qi > 5/3, when most of the 

duced through collisions. We investigate this parameter's produced cross-sectional area resides in the smaller grains 

effect in test runs exploring different values for s m i n ,coi (Fig.HSJ. 
(Fig. EJ). As can easily be seen, F max does not strongly 

vary with Smi , for Smin , col < 0.1 M m. There are two rea- 4 4 Qra]n chemjcal composition and j mpact strength 
sons for that: i) (3 decreases with decreasing sizes for grains 

smaller then ~ 0.1 (Fig.^), thus preventing them from The exact chemical composition of circumstellar disc ma- 

significantly contributing to the avalanche propagation; ii) terial is not well constrained and might in any case vary 

the broken power law for the debris size distribution, and from one system to the other. There is obs ervational ev- 

especially the flatter index qi — 1.5 for the smallest grains, idence for silicates, ices, and metals fe.g.. IPantin et all 
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Il997t iBouwman et al-lEoO^ . but their exact proportions 
in individual grains are difficult to estimate. Several de- 
tailed studies have addressed this issue for the specifi c 
/?Pic case (e.g.. lLi fc Greenber Jl998HPantin et al.lll997|) . 
but the estimates remain model dependent. Changes in 
grain compositions might affect the results in two ways: (i) 
compositional changes can lead to different values of /?, so 
that the grains experience different radiation pressure and 
as a consequence reach different outgoing (and impacting) 
velocities, and (ii) their collisional response properties can 
be significantly different. 

We first explore the role of grain porosities by varying 
this parameter between P = (compact grains, nominal 
case) and P = 0.8 (highly porous grains), with Qq remain- 
ing constant. This constant Qq prescription might seem 
counter-intuitive at first, since more porous grains should 
be expected to be more fragile, but it is, in fact, sup- 
ported by numerical experiments showing that porous tar- 
gets often prove more resistant than non-porous ones (e.g . , 
iFlvnn fc Durdi3l2004 iRvan et al.lfl99lt lLove et alJll99l , 
the reason being that impact shock waves are effectively 
dissipated by the pores. Figure ^2 shows that avalanche 
strength is maximum for the nominal case of compact 
grains (i^ ma x(p=o) — 210) and decreases for porous grains 

(-Pmax(P=0.8) — 70). 

We numerically explore the importance of chemical 
composition by performing runs for the 2 extreme cases 
of pure (compact) silicates (Mgo.gsFeo.osSiOa) and pure 
(compact) water ices. Here again, we take the possibly 
counter-intuitive constant Qq assumption, which is here 
again supported by experimental results showing that for 
target-projectile pairs of the same materia, !, ices can be as 
strong as silicates (e.g.. IRvan et al.lll999l) . Furthermore, 
compact ices and silicates of equivalent sizes have similar 
f3 values in the s > 0.1 /jm range (see Fig.^l. It is thus not 
surprising that our results show no significant difference 
between the pure-ice (-F ma x = 300) and pure-silicate runs 
(F max = 210). 

In a third set of runs we separately explore the Q* 
parameter, whose values for given grain compositions and 
dynamical conditions are still not well constrained by ex- 
periments or theoretical studies. The threshold energy es- 
timates for silicate-silicate and ice-ice collision s might vary 
betw e en ~ 10 6 and a few x 10 7 erg/g (e.g., Rvan et all 
Il999t iBenz fc Asphaudll999t iHolsapple et alJl2002|) . We 
explore Qq values between 10 6 and 10 8 erg/g and obtain 
strong variations in F max (Tig. ll2J) . For the lowest explored 
Qq value of 10 6 erg/g, we get F max — 10 4 , which is about 
50 times higher than in the nominal case (Qq = 10 7 erg/g). 

4.5. Field particle population 

As mentioned earlier, our reference field particle disc was 
assumed to be similar to the /3Pictoris system, for w hich 
the dust profile derived by lAugereau et al.l l)200ll) has 
been taken. Here we explore alternative profiles fFig. H3|) . 
Results show that F max does not strongly depend on the 



250 




J I I I I I I I L 



_L 



J I I L_ 



_L 



0.2 0.4 0.6 0.8 

porosity 

Fig. 11. Maximum amplification factor as a function of 
porosity for pure silicate grains. The value of the threshold 
energy, Q*, is assumed to be the same as in the nominal 
case. 
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Fig. 12. Maximum amplification factor vs. value of the 
threshold energy, Q*. Values for so = 1cm grains are de- 
noted on the axis. For the other sizes the threshold energy 
is given by Eq. 1111 



shape of the density distribution profile as long as the to- 
tal radial optical depth of the system (within 500 AU) rn 
remains the same. This result is in agreement with the 
simplified theory presented in Sect. [3 

On the other hand, we get drastic F max variations 
when changing the value of tii (regardless of the radial 
profile) . Figure ^] shows for example that increasing the 
number density by a factor of 5 leads to a value of 
-Fmax, which is a factor of ~ 1000 higher. This strong 
increase in .Fmax is in agreement with Eq. 1151 which 
predicts a strictly exponential growth with T|| , if (Np) 
is constant. However, in the simulations we find that 
(Np) weakly varies with tii through the empirical relation 
(Np) « 150 (T||/T|| !nom )~ - 45 . Plugging this expression for 
(Np) and (s(tj) 2 )'/(sl) ~ 1.5 into Eq.USwe get: 



ilinx — 6Xp 

= exp 



240 t,| (-3L 

5 - 3 (e=) 



0.55 



(16) 
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Fig. 13. Different test radial distributions for t_l (tii being 
constant). The thick soli d line is the dist r ibutio n for the 
nominal case, taken from lAugereau et al 1 !.-!()() ll). 
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Fig. 14. Maximum amplification factor as a function of 
tii . The open circles are the results of our simulations. 
The dashed line is the theoretical prediction fEa. 116ft . 

Thus tii proves to be the most efficient parameter for in- 
creasing Fmax by several orders of magnitude. This is a 
point of crucial importance when considering the absolute 
strength of an avalanche and its possible observability. 

5. Avalanche observability 

So far we have been concerned with the way an avalanche 
develops in a disc, in particular with how much dust can 
be created compared to the initial amount of released 
grains. This was quantified by the amplification factor 
F max - However, the crucial issue is under which condi- 
tions such an avalanche might become observable. In this 
respect, looking at F max is not enough. What matters here 
is the ratio between the luminosity of the avalanche par- 
ticles and that of the "field" particles, L av /L^. 

The value of L av /£d corresponding to the observabil- 
ity limit depends on several factors, such as the physical 
parameters of the system, observational conditions, and 
the observing devices' characteristics. Since our current 
study is not dedicated to a specific system, it is impossi- 
ble to give a precise criterion for avalanche observability. 



We shall thus adopt a simple and probably conservative 
criteria in which an avalanche is deemed observable when 
Lg, v I L& il 1 is reached at a given location in the disc. 

Most of the resolved debris disc images have been ob- 
tained in the visual or near-infrared (NIR) domains, domi- 
nated by scattered starlight. We shall thus focus here more 
specifically on scattered light luminosity. The amount of 
light scattered towards an observer coming from a given 
region of the disc is proportional to 



i(Ao) cx J J F*(\o,r)ns 2 



i (s,X )f(6)dn(s,r)dV, (17) 



where the integration is done over the spatial volume of 
the considered region and over the whole grain size range, 
s, with the grain number density n, the scattering coeffi- 
cient Qsca, and the scattering function f(9). F*(Ao,r) is 
the monochromatic star flux at the distance r from the 
star. 

For avalanche detection, the visible domain (~0.5 fim) 
is probably more favorable than the NIR (1 — 2 /jm). This 
is because avalanches consist mostly of submicron grains, 
which scatter very inefficiently at 1-2 /zm compared to big- 
ger grains a few microns in size (which is the average size 
for the "field" population). At the same time, in the vi- 
sual domain Q sca is nearly the same (within a factor of 2, 
depending on the exact chemical composition) for submi- 
cron and micron grains. Thus the ratio L av / is expected 
to be higher in the visual than in the NIR. For simplicity 
we assume that f(9) is only a function of the scattering 
angle 9 and that Q S ca is independent of the grain size. 
Although it is not exactly the case, this simplification can 
be considered as a reasonable starting point. 

Here we consider two extreme cases of disc orienta- 
tion, i.e., discs seen exactly edge-on and exactly pole- 
on. For the pole-on case, we determine that L av /Ld ~ 
T"_i_,av/T"_i_,field and thus use maps of the ratio between these 
vertical optical depths. For the edge-on case, we consider 
the synthetic midplane flux in scattered light, computed 
for different scattering functions. Since we do not know 
the exact optical properties of the circumstellar grains, 
two bracket cases have been considered for our calcula- 
tions: isotropic and forward scattering. For the forward 
scattering function we use an analytical approximation of 
the empirical / sca t for cometary d ust , obt ained from mea- 
surements of solar system comets l| Artvmowic dl 1997L and 
references therein): 



/scat(#) — fo 



0.3 



(0.2 + (9/2) 3 



1.4 



3.3 



0.2 



(18) 



In the following subsections, we investigate under 
which conditions avalanche-induced asymmetries might 
become observable for these two disc viewing angles. As 
appears from Figs. 1 1 51 1 81 these asymmetries consist of 
partial spiral or lumpy patterns in the face-on case, and of 
two-sided asymmetries, for which one side of the system 
becomes brighter than the other, in the edge-on case. We 
would like to point out that at other wavelengths, e.g., 
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Fig. 15. Top panel: Face-on case: Color-coded map of 
the ratio between the geometric surface densities of the 
avalanche grains and that of the "field" population for 
the nominal case. Bottom panels: Edge-on case: Midplane 
fluxes (arbitrary units) for the nominal case at avalanche 
maximum, edge-on orientation. The 2 solid lines indicates 
the total midplane fluxes ( "field" +avalanche) for each 
side of the disc (differences between the 2 sides are so 
small here that the 2 lines are almost indistinguishable). 
The dashed lines show the midplane fluxes for just the 
avalanche particles. Plot (a) corresponds to the forward 
scattering function and (b) to the isotropic case. 



infrared, the observability criteria (L av /Ld 1) might be 
reached for lower Ti.av/Ti, field ratios due to the fact that 
avalanche grains are expected to be hotter than field par- 
ticles, since their average size is about 10 times smaller 
then the average size of the initial disc population. 



Fig. 16. Same as Fig. 1151 but for collisionally weaker 
grains with Q^(sq = 1cm) = 10 6 erg/g (see Sect. I4.4f) . 



5.2. Larger amount of released dust Mo 

The most straightforward way of getting a more promi- 
nent avalanche is to increase the initially released amount 
of dust. As shown in Sect. 14.2.11 F max remains constant 
with varying Mq, so that the ratio L av /Ld increases lin- 
early with Mo. As a consequence, the release of as 10 22 g 
of dust would be required for the avalanche-induced lumi- 
nosities to become comparable to that of the rest of the 
disc. One might wonder however if a planetesimal shatter- 
ing releasing this large amount of dust is a common event 
(see discussion in Sect. I6.l)l . 



5.1. Nominal case 

As can be clearly seen in Fig. in the nominal case 
Lav/Ld never exceeds 10 -2 , neither in the edge-on nor in 
the head-on configuration. This value is far below our ob- 
servability criterion and the asymmetries induced by the 
corresponding avalanche would thus probably be unde- 
tectable by scattered light observations. 



5.3. Collisionally weaker grains 

As has been seen in Sect. 14.41 F max increases strongly for 
grains with lower specific energy values Q*. The lowest 
Q* value explored in Fig. EI Qo^icm) = 10 6er g/g> 
leads to rj_ av /Tl, field— 0.4-0.5. Thus, observability might 
be marginally reached when assuming the minimum shat- 
tering resistance for dust grains. 
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Fig. 17. Same as Fig. ^3 but for the case of a disc 5 times 
more massive than in the nominal case. 



5.4. Dustier systems 

The parameter exploration of Sect. 21 has clearly shown 
that the most efficient parameter for reaching high F max 
values is the field particles number density tii fEa. 
An obvious way of increasing rn is to assume a more 
massive disc, as has been done in Sect. 14.51 In terms of 
avalanche observability, we find that the observability cri- 
teria, L av /Ld ~ 1, is reached for a disc that is 4-5 times 
more dusty than in the nominal case. In this case, az- 
imuthal asymmetries become clearly visible in the face- 
on configuration and two-sided brightness asymmetries for 
the edge-on case fFig. H7[) . A massive dusty disc thus looks 
very promising from the point of view of avalanche obser- 
vation. 

Another way to reach higher values of tu is to keep 
the same total amount of dust, but distributed in a verti- 
cally thinner disc. In the nominal case the dusty disc has 
the characteristic width w (Eq. I1(J|) , a superexponential 
vertical profile (Eq. |5J), and a corresponding F max — 210. 



Assuming now a disc of thickness w* 
constant vertical profile 

dr = | r(R)/w* disc (R), if \z\ < w* disc 
dz \ 0, if \z\ > w^ isc , 



0.25w n 



with a 



(19) 
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Fig. 18. Same as Fig. 1151 but for a vertically thinner disc, 
described by Eq. ^| The total mass of the disc is the 
same as in the nominal case, but the asymmetries become 
prominent. 

we get F max = 2 x 10 4 for the same total amount of 
dust as in the nominal case. In this case, the number 
density of avalanche grains can even largely exceed that 
of the field particles (see Fig. I18|l. This density enhance- 
ment is azimuthally asymmetric due to the spiral struc- 
ture of an avalanche. If the disc is orientated face-on, then 
the azimuthal asymmetry persists for about 800 years. If 
the system is viewed edge-on then a two-sided asymme- 
try can be observed. Figure ED displays midplane fluxes 
( "field" -(-avalanche) for different azimuthal angles and 
scattering functions, clearly showing that, for favorable 
system orientations, one side gets significantly brighter be- 
cause of the avalanche. 



6. Discussions 

6.1. Probability of witnessing an avalanche event 

The numerical investigation of the previous sections has 
shown that collisional avalanches are a powerful and ef- 
ficient mechanism that naturally develops in debris discs 
after the breakup of a large planetesimal. However, in our 
nominal case of a f3 Pic-like system and Mq — 10 20 g of 
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dust initially released, the asymmetric features produced 
by the avalanche probably remain too weak to be observ- 
able in scattered light fSect. l5"T|) . This result should, how- 
ever, be taken with great care since our parameter explo- 
ration has shown that avalanche strength strongly depends 
on several critical and often poorly constrained parame- 
ters. The first set of parameters is those linked to the 
initial breakup event. Here we obtain the intuitive result 
that higher amounts of initially released dust leads to more 
powerful avalanches (see Sect. 15.2(1 . with the avalanche 
strength scaling linearly with M n . This is not uniq ue to 
the avalanche mechanism: IKenvon fc BromlevI l|2005h find 
a similar dependence when only considering the signature 
of the cloud of primary debris produced immediately after 
the planetesimal breakup. What distinguishes our results 
from studies in which only dust released at impact is con- 
sidered is that avalanches strongly depend on the number 
density of dust in the disc. Section l4~51 has indeed shown 
that the global optical depth of the dust disc rn is the pa- 
rameter avalanche development depends most on, the de- 
pendence being close to an exponential. We have seen that 
other parameters, mostly related to the way the physical 
response of grains to collisions is modeled, might also lead 
to observable events when stretched to the extreme values 
that were numerically explored here. This is in particu- 
lar the case for Q*, for which very low ~ 10 6 erg/g values 
might lead to powerful avalanches. 

We shall however leave these "technical" parameters 
aside to focus on the 2 parameters directly related to 
the system's properties themselves, i.e., the optical depth, 
both t± and T|| , and the initial amount of dust released 
M , and derive an order-of-magnitude estimate for the 
probability of witnessing avalanche events as a function of 
these parameters. From the results of Sect.0] the L av /Ld 
criterion for observability might be written 



F 



max(nom) 



Mo 

10 20 g 



is 100, 



(20) 



which is equivalent to saying that the luminosity ratio 
between avalanche and field grains should be at least 100 
times higher than in the nominal case (for which L av /Ld ~ 
10~ 2 ). Section l4.2.1l has shown that P ma x is independent of 
Mo, so that in our approximation F max is only a function 
of Ti| , and this Ty dependence is given by Eq. 1161 Thus, 
Eq. |201 reduces to 



exp 



5.3 



| ,nom 



0.55 



M 

10 20 g 



TL,nom 



Z 2 x 10 4 



(21) 



which gives a direct link between a given disc density 
(tj_ and T||) and the minimum mass of released dust 
able to produce a visible avalanche in such a disc (the 
denser the disc, the smaller the corresponding Mo value). 
The other important issue affecting witnessing probabili- 
ties is the duration of an avalanche. Our simulations show 
that the typical lifetime of an avalanche-induced pattern 



is t av ~ 10 3 yrs. With this value and Eq. |^ one can es- 
timate the probability P b s of witnessing an observable 
avalanche event in a given disc: 



-fobs — 



£shatt(M ,Tx) 



(22) 



where t s hatt(M ,T) is the average time between 2 shatter- 
ings producing Mq of dust in a disc of average optical 
depth tj_, with Mo, n and t± satisfying Eq. [21 As sug- 
gested in Sect. 13.51 we consider that the object releasing 
M of dust ha s a mass Mpr c± lOMn . From unpublished 
results of the iThebault et al.l l|2003l) simulations of col- 
lisional rates and outcomes in the inner /3Pic disc, we 
determine that the approximate timescale for the shatter- 
ing of a Mpb = IOMq object to occur in the innermost 
< 50 AU (the typical location for the initial shattering 
events considered in our simulations) of a /3Pic like sys- 
tem is i sh att ^ 150[(10Af )/10 21 g] 1 - 25 yrs. Since, for sys- 
tems with similar spatial distributions, the frequency of 
collisional events is proportional to the square of a sys- 
tem's total mass, we get the empirical relation: 



shatt(Mo,Tj_) 



150 



( M 



V 10 20 i 



1.25 



yrs, (23) 



where we implicitly assume that the system's spatial dis- 
tribution is the same as in the nominal case, so that the ra- 
tio between two systems' total masses is equal to the ratio 
ti/ti nom anywhere in the disc. This equation should of 
course be regarded as giving a very rough estimate, since 
*shatt(Mo,Tj_) depends on many poorly constrained param- 
eters, like the number density of planetesimals and their 
average kinetic energy at impact. Equation 1231 however, 
gives the global trend of the way t s hatt(M ,Tjj increases 
with Mq. Taking the lowest Mo value satisfying Eq. |^ 
and plugging it into Eq. [231 we get, from Eq. 1221 



Pobs «3xl0~ 5 



0.75 



exp 6.6 



0.55N 



.(24) 



Equation 1241 indicates that P ODS — 0.03 for the nominal 
case field particle disc, which means that we have about 
a 3% chance of witnessing the avalanche caused by the 
breakup of a Mpb = 10M -10 23 g object (M = 10 22 g 
being the smallest released dust mass able to trigger a vis- 
ible avalanche for such a disc, as given by Eq. 121 [I . This 
makes it a rather unlikely event, although it cannot be 
completely ruled out. Nevertheless, slightly denser discs 
(i.e., higher t_l and rii) can easily raise P bs up to 1. As 
a matter of fact, the dependence on rn is so sharp that 



P 



obs 



1 is obtained for tii ~ 2.1t\ 



|| ,nom 



0.046. We thus 



see that a (3 Pic-like system is below, but not too far from 
the limit for which chances of witnessing an avalanche 
are high, especially when considering the uncertainties re- 
garding avalanche strength due to its dependence on sev- 
eral poorly constrained parameters related to the collision- 
outcome prescription (also keeping in mind that higher T|| 
values could alternatively be achieved for a thinner disc 
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of the same dust mass (see Sect. 15.40 ). Moreover, our 
L&v/Ld > 1 criterion for observability is probably too 
conservative, and avalanche-induced patterns might be de- 
tectable for lower luminosity excess values. Taking, for ex- 
ample, L a v/id > 0.1 would raise the detection probability 
to ~ 45% for a (3 Pic-like system. 

6.2. Avalanches in observed systems, perspectives 

We defer a detailed application of our model to specific 
circumstellar discs to a future study. However, the present 
results can already give a good idea of the typical profile 
for a "good" avalanche-system candidate 

Our numerical exploration has shown that structures 
that are the most likely to be associated with avalanche- 
events have two-sided asymmetry for discs viewed edge-on 
and open spiral patterns for discs viewed pole-on or at in- 
termediate inclinations. An additional requirement is that 
these discs should be dust-rich systems, with a dustiness 
at least equal to, and preferably higher than that of /3-Pic. 
Note also that our model makes an additional prediction, 
i.e., that avalanche affected regions should consist of grains 
significantly smaller than the "field" particles in the rest 
of the disc. If the blow-out radius of grains is of the order 
of the wavelength of the observed light, then this should 
translate into color differences between avalanche (bluer) 
and non-avalanche (redder) regions. 

In this respect, one good edge-on candidate might 
be the recently discovered HD 32297 system, which ex- 
hibits a strong two-s ided asymmetry. As reported by 
ISchneider etafl j2005|) and iKalasI (2005), this system is 
a (3 Pic analog with its SW sid e significantly brighter than 
the NE one within ~ 100 AU Schneider et afll2005|) and 
possibly outside 500 AU ftCalad 12005^ Such a two-sided 
asymmetry would be compatible with the ones obtained in 
our simulations (as shown f or example in the bottom panel 
of Fig. UHJ). Furthermore, iKalad (|2005^ also reported a 
color asymmetry between the two sides, with the brighter 
one (SW) being significantly bluer. This seems to indicate 
that this side is made of smaller, possibly submicron grains 
l)KalasH2oT)5T) . As previously discussed, this is what should 
be expected for an avalanche-affected region. However, 
an alternative scenario, like the colli sion with a clump 
of interstellar medium proposed by jKalas 2005), might 
also explain the HD 32287 disc structure. Future imaging 
and spectroscopic observations are probably needed before 
reaching any definitive conclusions. 

Among all head-on observed systems, the one display- 
ing the mo st avalanche- like str ucture is without doubt 
HD 141569 ( Clampin et al.l20'03T) . with its pronounced spi- 
ral pattern. Furthermore, the disc's mass, significantly 
higher than /3-Pictoris, makes it a perfect candidate in 
terms of witnessing probabilities. Of course, avalanche is 
not the only possible scenario here, and several alterna- 
tive explanations, like an eccentric bound planet or stel- 
lar compani on, or a stellar flyby have already been pro- 
posed (e.g., lAueereau fc Papaloizou] 120041: IWvattl 120051 



lArdila et al.ll2005h . One should, however, be aware that 
this system is stri ctly speaking no t a "s tandard" debris 
disc as defined by lLagrange et alJ l)2000h and as consid- 
ered in the present study. Indeed, several studies seem to 
suggest the presence of large amounts of primordial gas 
ijZuckerman et alJll995l lArdila et al.ll2005() . 

Gas drag effects have been left out of the present study 
on purpose, mainly because, in the strict sense of the term, 
debris discs are systems where dust dynamics is not dom- 
inated by gas friction (|Lagrange et alJEoOfl) . Moreover, 
the correct description of dust-gas coupling adds several 
additional free parameters (gas density and temperature 
distributions, etc.) and requires a full 2-D or 3-D treat- 
ment of gas by far exceeding the scope of the present pa- 
per. However, the issue of avalanches in a gaseous medium 
might be a crucial one for those systems that are most fa- 
vorable for avalanches, i.e., discs more dusty than /3-Pic, a 
system which is already at the upper end of de bris-discs in 
terms of dustiness (e.g.. ISnangler et al"ll200lj) . Such more 
massive systems should fall into a loosely defined category 
of "transition" discs between T-Tauri or Herbig Ae proto- 
planetary syste ms and "proper" de bris discs (see for ex- 
ample Sect. 4 of lDutrev et alJ l2004). For such systems (of 
which HD141569 is a typical example), which are younger 
than the more evolved debris discs, risks (or chances) of 
encountering large amounts of remaining gas are high. 
This crucial issue will be the subject of a forthcoming 
paper. 

7. Summary 

This paper presents the first quantitative study of the col- 
lisional avalanche process in debris discs, i.e., the chain 
reaction of dust grain collisions triggered by the initial 
breakup of a planetesimal-like body. We have developed a 
code that allows us to simultaneously follow both spatial 
and size distributions of the dust grain population, which 
is collisionally evolving because of impacts caused by small 
particles blown out of the system by the star's radiation 
pressure. Our results can be summarized as follows: 

1. Collisional avalanches propagate outwards leaving 
a characteristic spiral-shaped pattern in the system. 
Depending on the system's orientation, these patterns 
might appear as open spirals or lumpy structures (face- 
on geometry) or a two-sided brightness asymmetry 
(edge-on case). In a j3 Pic-like disc, an avalanche lasts 
for about 10 3 years. 

2. The strength of an avalanche depends linearly on the 
mass of the initially shattered object, but nearly ex- 
ponentially on the optical depth of the dust disc in 
which it propagates. The disc's dustiness is by far the 
most crucial parameter here, making dusty discs much 
more favorable cases for avalanche propagation than 
tenuous ones. 

3. We define a conservative criterion for avalanche observ- 
ability, from which we infer a relation between a given 
disc density and the minimum mass of the object that 
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has to be shattered to reach observability. When cou- 
pling this relation to estimates for catastrophic disrup- 
tion probabilities among planetesimal-objects in debris 
discs, we are able to derive a first-order estimate for 
the probability of witnessing an observable avalanche 
event in a given debris disc. For our reference (3 Pic-like 
system it is of a few percents, but probabilities rapidly 
increase for slightly denser systems. 
4. Modeling of dustier young transitional discs may re- 
quire the inclusion of gas drag, which may change both 
the morphology and the strength of the avalanche. 

Appendix A: Superparticles' structure 

All SPs are cylinders with variable height h sp and constant 
radius r sp . Their geometrical centers always stay in the 
disc midplane. For most of the runs we take r sp = 5 AU. 
Test runs showed that an avalanche development weakly 
depends on the value for r sp . The maximum amplification 
factor for runs with r sp in the 4 — 8 AU range differs by less 
then ~ 10%. SPs are modeled in different ways depending 
on the physical origin of the grains they represent. We 
distinguish between 3 types of SP. 



A.l. Field SP 

All SPs that represent the initial structure of the dusty 
disc (i.e., non-avalanche SPs) have a superexponential ver- 
tical density profile 



sp,i f-0.5h SPii f( Z )d z 



f(z), 



(A.l) 



where A^ gr .i is the total number of grains in a given SP, 
r sp andh sp are the SP's radius and height, and f(z) is the 
SP profile function, which reads 

'<"-- » H^mT)- (A ' 2) 

where p z = 0.7 and w sPl i is the SP's width, which depends 
on the distance to the star, R, and is equal to the disc 
width (Eq. HUJ for field SPs. The height of a field SP is 



(R) = 2Qw SPj i(R). 



(A.3) 



All field SPs standing for the biggest grains (1cm) 
have circular orbits. The SPs radial distribution and the 
number of grains in each SP is chosen in accordance with 
the be st-fit parent body distribution from lAueereau et alJ 
( 2001). All other grains are assumed to be produced from 
these biggest grains following the power law size-frequency 
distribution of Eq. [S] The number of SPs in each size bin 
is taken such that in the steady-state configuration there 
are 2-5 overlapping SPs of the same grain size at any given 
location in the system. 

A. 2. Initially released planetesimal debris 



(1998| show that the fragment velocities depend on grain 
sizes and that for small grains they are about a few per- 
cent of the impact velocity, with Vf r w 0.01 — 0.1-Ui mpac t. 
Fragments produced from a body on a Keplerian orbit 
would spread in the vertical direction, and the height of 
the layer can be estimated as 



h 



sp 



^impact 



-R. 



(A.4) 



For most of the runs we take ffr/fimpact = 0.1. 

If the debris have velocities with isotropic distributions 
around the center of mass then the vertical distribution of 
grains of the same f3 can be approximated as a constant. 
Thus these SPs are assumed to have no internal density 
structure and constant grain density. 



gr ' 1 ~ — 2 A 



(A.5) 



' sp,i'"sp,J 



which corresponds to f(z) = 1 (i.e., io sp ,i = oo in Eq. 

rot . 

A.3. SP created due to collisions 

As mentioned in Sect. [3| when two SPs are passing 
through each other, a fraction of their grains can be de- 
stroyed producing new and smaller grains, which are then 
combined into new SPs. The number of grains that are de- 
stroyed in each SP, if the relative velocity is high enough 
to reach catastrophic fragmentation (see Sect. l3.4|l . is cal- 
culated as 



7V„ r ,_ = 



(A.6) 



where z = 0.5 min(/i sp ,j, /i S p,j)> a = 7r(s grii + s gr j) 2 /4 
is the collisional cross-Sect, for grains with physical sizes 
Sgr(i,j)) Vrd is the relative velocity of the grains, A ovcr is 
the overlapping area of the two SPs viewed top-on, and t co \ 
is the time while the SPs are passing through each other. 
The time dependences in Eq. I A. 61 might be neglected if 
a reasonably small "collisional" time step is taken for the 
calculations. We adopt At co \ = 0.02 (in units of the orbital 
period at 20 AU) for most of the runs. Simulations with 
much smaller collisional time steps (e.g., At co \ = 0.005) 
do not lead to a significant improvement of the results. 
The debris velocities and the size-frequency distribution 
for the newly created grains are identical to the values 
that one would get, considering a collision between N gIt - 
pairs of grains with sizes s gr ,(ij) and relative velocity v rc \ 
(Sect. EH . 

The structure of the new SPs (one SP for each size bin), 
which are produced after a collision, is obtained through 
the following equations. The initial height of the SP is 
equal to 



When a planetesimal is shattered the fr agments are ere 
ated with an initial spread in velocities. iRvan fc Meioshl f k (z) 



h sp ,k = min(h sPli , h spj ). (A. 7) 

The SP vertical profile fk(z) is calculated as 

fi(z)fj(z), (A.8) 
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which corresponds to a SP of width (Eq. IA.2fl 



the probability proportional to 



iV 2 



Choosing posi- 



W S p,k 



1 -1/P* 



(A.9) 



In the case of dust production by SPs with super- 
exponential profiles (Eq. IA.2|) . the new SPs will have a 
smaller width than the two "colliding" SPs. Keeping the 
same width would allow us to overlook the fact that "real" 
grains in the SPs have vertical velocities and so will the 
produced debris. This component can be neglected if we 
only consider collisional outcomes, but it should be taken 
into account for calculations of the density structures. In 
our calculations this effect is taken into account by increas- 
ing the newly-created-SPs' width and height with time. 
The increase rate is taken to be 



At,, 



(A.10) 



where At w is equal to \ of the orbital period at the lo- 
cation where the SP is created. The same kind of time 
dependence is applied to the h sp growth, namely 



-h so 



At h 







if ^sp ^ ^max) 



(A.ll) 



where h max — ZioWdisc, th = h$t w , and ho is a constant 
parameter. The results for ho = 0.3 and ho = 3 differ by a 
factor of 2, which is acceptable for our order-of-magnitudc 
calculations. 

Appendix B: Recombining SPs 

To keep the total number of SPs manageable (we can trace 
the evolution of about one million of them) we have devel- 
oped an algorithm that allows us to recombine SPs with 
similar parameters (grain size, velocities, positions in the 
disc). This allows us to avoid a too fast increase of the to- 
tal number of SPs while keeping it to a value large enough 
from a statistical point of view. 

The merging procedure is applied only between SPs 
standing for grains of the same size. The proximity con- 
dition is obtained by dividing the disc plane into a 2-D 
space grid with the cell size equal to the SP's diameter. 
For each cell we list all SPs whose centers are located 
within the cell. Thus a given list contains SPs from the 
same spatial volume and with the same grain size. For 
each SP from a given list we calculate the "normalized" 
velocity v* — v/v\, cp (R), where v is the SP velocity and 
R is the distance from the star to the SP center. If sev- 
eral SPs fall into one velocity bin (of width dv* = 10~ 2 ), 
they are combined into one SPs with the number of grains 
being equal to the sum of N SI of the combined SPs. The 
vertical structure of the new SP is obtained through the 

averaged values w sv = „ and h sv = £ ■ 

The velocity and the position of the new SP are chosen to 
be equal to those of one SP, randomly chosen from the list 
of the recombined SPs. This SP is chosen randomly with 



tions and velocities this way has the advantage of not 
introducing new trajectories for SPs, and this makes us 
more certain about the general shape of an avalanche. We 
have also tested another procedure for SPs recombination 
in which all SPs from the same velocity bin are recombined 
into one SP with position and velocity such that the total 
angular momentum and kinetic energy are preserved. In 
both cases we obtain similar results. 

The recombination described above makes a significant 
reduction of the total number of avalanche SP spossible. 
However, an additional optimization procedure has been 
implemented to speed up the calculations. The idea is the 
following: The dust production rate is approximately pro- 
portional to the number density of avalanche grains (mul- 
tiplied by the number density of the field population) . As a 
consequence, regions with the lowest density of avalanche 
grains (of a given size) cannot make a significant contri- 
bution to avalanche dust production. Thus these regions 
are not very interesting for our simulations and there is no 
need to do very accurate representations of the avalanche 
grains' population there. A very rough representation is 
enough here. The question is to determine which number 
density values should be considered as "unproductive". 
Since midplane densities for the field population vary by 
a factor ~ 100 throughout the disc, regions where the den- 
sity of avalanche grains (of a given size bin) are less than 
1/1000 that of the most dense regions will be considered 
as " not important" . In these regions all SPs from the same 
space cell (the proximity condition as defined above) rep- 
resenting grains of the same size bin are combined into 
one SP regardless of their velocities. Test runs have proved 
the efficiency of this optimization procedure: i) the gen- 
eral shape of an avalanche and the amplification factor 
evolution are not modified by this procedure, since the 
"important" regions are not affected; ii) the total num- 
ber of SPs is significantly reduced, significantly speeding 
up the calculations. However, this procedure is limited to 
cases for which there is a significant density gradient. At 
later stages of avalanche propagation (^ 12 orbital peri- 
ods) it is inefficient. Figure lB~D shows the evolution of the 
total number of SPs with time. 
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